﻿using System;
using System.Collections.Generic;
using OSGeo.OSR;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using System.Windows;

namespace GPVC
{
    class GetPixelValuesByCoordinates
    {
        private static List<RasterIO.GeoTiffInfo> geoTiffInfoList = new List<RasterIO.GeoTiffInfo>();
        public static List<RasterIO.GeoTiffInfo> GeoTiffInfoList { set { geoTiffInfoList = value; } }

        private static ExcelIO.ExcelInfo excelInfo = new ExcelIO.ExcelInfo();
        public static ExcelIO.ExcelInfo ExcelInfo { set { excelInfo = value; } }


        private static Tuple<double, double> LonLat2XY(SpatialReference gcs, SpatialReference pcs, double[] lonlat)
        {
            CoordinateTransformation ct = new CoordinateTransformation(gcs, pcs);
            Tuple<double, double> xy = new Tuple<double, double>(lonlat[0], lonlat[1]);
            //
            try
            {
                ct.TransformPoint(lonlat);
                //
                xy = new Tuple<double, double>(lonlat[0], lonlat[1]);
            }
            catch (Exception errInfo)
            {
                var dlgResult = MessageBox.Show(errInfo.Message,
                                 "错误", MessageBoxButton.OK, MessageBoxImage.Error);
                if (dlgResult == MessageBoxResult.OK) return xy;
            }
            //
            return xy;
        }

        private static Tuple<double, double> XY2LonLat(SpatialReference gcs, SpatialReference pcs, double[] xy)
        {
            CoordinateTransformation ct = new CoordinateTransformation(pcs, gcs);
            Tuple<double, double> lonlat = new Tuple<double, double>(xy[0], xy[1]);
            //
            try
            {
                ct.TransformPoint(xy);
                lonlat = new Tuple<double, double>(xy[0], xy[1]);
            }
            catch (Exception errInfo)
            {
                var dlgResult = MessageBox.Show(errInfo.Message,
                                 "错误", MessageBoxButton.OK, MessageBoxImage.Error);
                if (dlgResult ==  MessageBoxResult.OK) return lonlat;
            }
            //
            return lonlat;
        }

        private static Tuple<double, double> XY2RowCol(double[] geoTransform, double[] xy)
        {
            Matrix<double> matrixA = DenseMatrix.OfArray(new double[,] { { geoTransform[1], geoTransform[2] },
                                                                         { geoTransform[4], geoTransform[5] } });

            DenseVector vectorB = new DenseVector(new double[] { xy[0] - geoTransform[0], xy[1] - geoTransform[3] });

            Vector<double> rowcolVector = matrixA.LU().Solve(vectorB);
            double row = rowcolVector[1];
            double col = rowcolVector[0];
            //
            Tuple<double, double> rowcol = new Tuple<double, double>(row, col);
            //
            return rowcol;
        }

        private static Tuple<double, double, double> RowCol2XY(double[] geoTransform, int row, int col)
        {
            double x = geoTransform[0] + row * geoTransform[1] + col * geoTransform[2];
            double y = geoTransform[3] + row * geoTransform[4] + col * geoTransform[5];
            //
            Tuple<double, double, double> xyz = new Tuple<double, double, double>(x, y, -1);
            //
            return xyz;
        }

        public static double GetPixelValues(double[] bandArr,
                                            int width,
                                            bool hasPcs,
                                            OSGeo.OSR.SpatialReference gcs,
                                            OSGeo.OSR.SpatialReference pcs,
                                            double[] geoTransform,
                                            double[] coordinates,
                                            string coordinatesType)
        {
            double rasterValue = double.NaN;
            int row = -1;
            int col = -1;
            //
            if (coordinatesType == "ROWCOL")
            {
                row = Convert.ToInt32(coordinates[0]);
                col = Convert.ToInt32(coordinates[1]);
            }
            if (coordinatesType == "LONLAT" && hasPcs == false)
            {
                double[] srcCoordinates = { coordinates[0], coordinates[1], 0 };
                Tuple<double, double> rowcol = XY2RowCol(geoTransform, srcCoordinates);
                row = (int)rowcol.Item1;
                col = (int)rowcol.Item2;
            }
            if (coordinatesType == "LONLAT" && hasPcs == true)
            {
                double[] srcCoordinates = { coordinates[0], coordinates[1], 0 };
                Tuple<double, double> tmpXy = LonLat2XY(gcs, pcs, srcCoordinates);
                double[] xy = { tmpXy.Item1, tmpXy.Item2 };
                if (xy[0] != srcCoordinates[0] && xy[1] != srcCoordinates[1])
                {
                    Tuple<double, double> rowcol = XY2RowCol(geoTransform, xy);
                    row = (int)rowcol.Item1;
                    col = (int)rowcol.Item2;
                }
            }
            if (coordinatesType == "XY")
            {
                double[] srcCoordinates = { coordinates[0], coordinates[1], 0 };
                Tuple<double, double> rowcol = XY2RowCol(geoTransform, srcCoordinates);
                row = (int)rowcol.Item1;
                col = (int)rowcol.Item2;
            }
            if (row < 0 || col < 0) return rasterValue;
            else rasterValue = Convert.ToDouble(bandArr.GetValue(row * width + col));
            //
            return rasterValue;
        }


    }
}
